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ABSTRACT 


This  thesis  considers  two  models  for  the  computation  of 
position  finding  confidence,  one  of  which  utilizes  a  bivariate 
normal  distribution  of  region,  and  the  other  a  chi-square 
distribution.   The  two  models  are  based  on  different 
assumptions,  these  are  explained  and  explored.   A  computer 
simulation  model  is  presented  which  utilizes  both  position 
finding  models  under  varying  conditions. 
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I.   INTRODUCTION 

Bearings  are  taken  on  an  emitter  (target)  of  unknown 
geographic  location  by  n  stations  of  known  location.   The 
best  point  estimate  (BPE)  of  the  target  is  then  calculated, 
along  with  a  confidence,  or  search,  region  within  which  the 
target  may  be  said  to  be  located  with  a  given  probability. 
This  constitutes  the  classic  radio  direction  finding  position 
location  problem. 

This  classic  problem  is  part  of  a  general  class  of 
position  estimation  problems  whose  basic  approach  is  identical, 
which  is  the  calculation  of  a  point  location  from  a  number  of 
lines  of  position. 

Variant  problems  in  this  class  occur  when  lines  of  position 
are  established  from  ranges,  or  relative  times  of  arrival  of 
an  emission,  or  when  the  object  to  be  located  takes  the 
measurements  on  the  reference  stations. 

While  this  paper  examines  the  classic  problem,  much  of  the 
methodology  applies  to  this  more  general  class  of  problems. 

The  BPE  is  the  primary  result  of  a  position  location 
solution.   The  confidence  region  amplifies  this  information 
in  two  ways.   The  size  of  the  confidence  region  serves  as  a 
measure  of  the  reliability  of  the  BPE,  a  larger  region 
implying  a  more  tenuous  location.   The  region  also  describes 
an  area  in  which  the  target  will  be  located  with  specified 
confidence.   The  geometric  details  of  this  region  provide  the 
more  likely  locations  for  search. 


II.   COMPUTATION  OF  CONFIDENCE  REGIONS 

A  typical  [6]  sequence  of  operations  for  a  position 
estimation  by  computer  is  as  follows:   A  first  point  esti- 
mate is  made  and  a  plane  surface  tangent  to  the  earth  is 
established  at  that  point.   Remaining  computations  are  made 
in  the  plane  rather  than  the  surface  of  the  earth.   Unusable, 
or  wild,  bearings  are  rejected  and  a  BPE  made  utilizing  the 
remaining  bearings.   A  confidence  region  is  computed  using 
the  same  bearings. 

More  than  one  mathematical  method  is  available  for  most 
of  the  steps  [3,  4,  8,  9,  11] .   The  two  methods  used  herein 
for  confidence  region  computation  are  designated  the 
bivariate  normal  (BVN)  and  the  chi-square  (x  ) • 

The  Bivariate  normal  region  [3,  4]  is  the  most  commonly 
used.   This  region  is  the  family  of  concentric  ellipses 
ax2  -  2bxy  +  cy2  =  -2  log  (1  -  P) 

where  P  =  probability  that  the  target  is  found  within  the 

ellipse 

x,  y  =  coordinates  in  local  (to  the  target)  reference 
system 
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0 •   =  bearing  observed  by  station  j 

9         •  th   •  • 

a*   =  variance  of  displacement  of  ]    line  of  bearing 

n    =  number  of  bearings 

The  bearings  from  the  stations  are  assumed  to  be  normally 
distributed  with  a  mean  of  zero,  and  are  displaced  parallel 
to  the  true  bearing  at  the  site  of  the  BPE.   The  confidence 
region  lies  in  a  plane  tangent  to  the  surface  of  the  earth 
at  the  BPE  as  calculated  by  the  least  squares  method  [3,  4]. 

The  BVN  confidence  region  requires  several  assumptions: 

(1)  The  position  lines  are  great  circles  on  the 
surface  of  the  earth. 

(2)  Measurement  error  is  Gaussian  with  mean  zero. 

(3)  Displacement  of  position  lines  from  the  true 
position  line  is  parallel  to  the  true  position 
line. 

(4)  The  earth  is  flat  in  the  vicinity  of  the 
confidence  region. 

These  last  two  assumptions  represent  approximations  that 
are  made  to  simplify  the  calculation  of  the  confidence  region. 
Thus,  the  calculated  region  is  an  approximation.   The  parallel 
bearing  displacement  is  an  approximation  of  the  actual  dis- 
placement of  a  point  on  the  bearing  by  a  bearing  error.   From 
this,  the  probability  distribution  of  the  points  is  developed 
[4] .   Placing  the  resultant  curves  in  a  plane  adds  an 
additional  distortion,  which  grows  with  the  size  of  the 
confidence  region  [7]. 
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The  chi-square  confidence  region  [3,  10]  is  the  family 
of  concentric  curves 


n  9 


x*  =   e  -4-.  (e .  -  6.) 


v;here 

S .   =  bearing  variance  of  station  j 
3 

0 .   =  bearing  observed  by  station  j 

3    =  true  bearing  from  station  j  to  a  point 
j 

(x,y)  in  the  curve 
n    =  number  of  bearings 
X   has  a  chi  square  distribution  with  n  degrees  of 
freedom.   Of  the  above  listed  assumptions  required  by  the 
BVN  method,  only  (1)  and  (2)  are  required  for  the  x   region. 

The  region  lies  on  the  surface  of  the  earth  centered  at  the 

2 

minimum  X   point,  rather  than  in  a  plane. 

Both  methods  require  the  estimated  variance  of  the 
bearings  to  be  used  in  the  calculation.   The  effect  of  an 
incorrect  estimation  could  be  much  larger  than  the  errors 
of  approximation.   Estimates  in  practice  are  a  function  of  a 
stations  past  performance .plus  the  direction  finder  operator's 
evaluation  of  the  reliability  of  each  bearing  he  observes  and 
reports . 

The  assumption  of  Gaussian  bearing  distribution  eliminates 
the  case  of  wild,  or  outlying  bearings  which  occur  in  practice 
largely  through  operator  error.   As  a  Gaussian  distribution 


2 
is  basic  to  both  the  BVN  and  x   methods,  elimination  of  wild 

bearings  from  the  computation  is  necessary,  though  the  theory 

of  each  method  would  be  affected  by  their  assumed  existence. 

Wild  bearing  rejection  methods  are  available  [6,  7],   The  use 

of  such  methods  is  usually  limited  to  four  or  more  bearing 

situations.   With  two  bearings,  no  basis  for  detection  of 

wild  bearings  exists.   For  most  three  bearing  situations,  when 

the  existence  of  at  least  one  wild  bearing  is  indicated,  it 

is  impossible  to  identify  it,  or  them.   With  four  bearings, 

one  wild  bearing  can  be  seen  to  miss  the  BPE  area  formed  by 

the  other  three. 

In  practice  [6,  7]  the  Outlier  detection  and  rejection 
method  is  used,  but  a  Gaussian  distribution  is  still  assumed 
for  those  remaining. 

Daniels  and  Vajda  [3]  object  to  the  use  of  x  regions  on 
the  grounds  that  it  may  occur  that  no  points  can  be  found  to 
satisfy  the  resultant  equations. 

Read  [10]  examines  the  simplifications  required  in 
calculating  the  BVN  regions,  and  compares  the  two  methods. 
The  x   regions  are  expected  to  be  too  small,  and  the  BVN  too 
large,  though  the  full  effect  of  the  simplifications  by  Kukes 
and  Starik  [4]  is  not  known. 

Comparison  of  the  two  methods  by  computer  plotting 
techniques  shows  considerable  difference  in  size.   Additionally, 
the  effect  of  the  approximations  used  in  the  BVN  method  is 

apparent.   While  BVN  region  is  elliptical  in  all  cases,  the 

2 
X   region  may  assume  an  egg  shape. 


A  barrier  to  the  use  of  x   regions  is  the  lack  of  a  way 
of  reporting  them.   A  BVN  region  can  be  described  by  its 
center,  angular  orientation,  and  length  of  semi  axes.   No 
such  simplifying  parameters  exist  for  x   regions.   One 
possible  approach  would  be  to  attempt  a  conversion  into  a 
partially  equivalent  rectangle,  square,  or  circle  as  is  often 
done  for  BVN  regions  [1,  6,  7].   Since  these  regions  are 
always  larger  than  the  computed  region,  contain  areas  of 
lesser  probability,  and  omit  areas  of  greater  probability, 
any  advantage  of  using  one  method  over  the  other  would 
probably  be  lost  thereby. 

A  thorough  investigation  of  the  two  methods  is  required, 
in  order  to  determine  the  actual  probabilities  represented 
by  each  type  of  confidence  region.   The  differences  indicate 
that  at  least  one  of  the  methods  has  error.   A  computer 
simulation  was  written  to  compare  them. 


III.   COMPUTER  SIMULATION  MODEL 

The  appended  computer  program  provides  a  means  for 
simulation  of  the  two  confidence  region  methods  for  identical 
input  data. 

Inputs  to  the  program  are  station  latitude,  longitude, 
and  standard  deviation  of  bearing  distribution;  target 
latitude  'and  longitude. 

The  true  bearings  and  distance  for  each  station  and 
target  are  computed  and  a  bearing  error  generated  and  added. 
This  portion  contains  parts  of  the  Bergh  program  in  Pope  [8] 
and  is  a  program  used  by  LT  Paul  Winkler,  modified  here  to 
include  more  general  positions. 

A  first-point  estimate  of  the  position  is  made  by  Pope's 
[9]  vector  method,  and  the  final  fix  by  the  least  squares 
method  of  Kukes  and  Starik  [4]  in  a  plane  tangent  to  the 
earth  at  the  first-point  estimate.   This  portion  was  used  by 
Lunde  [5]  but  contained  errors  since  corrected. 

A  bivariate  normal  confidence  region  is  computed  by  the 
method  of  Kukes  and  Starik,  and  orientation  angle  and 
length  of  major  and  minor  simi  axes  generated.   Additionally, 
output  by  graphic  plot  is  available. 

A  x   confidence  region  is  then  computed  and  a  graphic 

2 
plot  generated.   No  parameteric  output  for  the  x   region  is 

available,  as  no  method  yet  exists  for  describing  the  shape 

of  the  region.   This  portion  was  originally  written  by  LT 

Winkler,  and  since  modified  to  include  more  general  positions 
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The  subroutines  GAMA,  for  the  normalized  incomplete  gamma 
function,  and  MTRMAP ,  the  graphic  plot,  are  from  the  library 
of  the  W.  R.  Church  Computer  Center  at  the  Naval  Postgraduate 
School . 

The  graphic  plots  are  on  a  grid  measured  in  degrees  of  a 
great  circle,  resembling  a  Mercator  projection  at  the  equator 
Hence,  visual  distortion  of  the  regions  occurs  as  their 
location  approaches  the  poles,  but  the  geographical  points 
of  the  x   regions  are  correct.   The  geographical  points  of 
the  BVN  regions  are  correct,  except  for  the  distortion  caused 
by  the  flat  earth  assumption. 
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IV.   SUGGESTED  SIMULATION  PROCEDURE 

An  undetected  error  in  previous  work  invalidated  all 
data  generated.   Time  constraints  prevented  regeneration  of 
the  data  after  correction  of  the  error. 

The  behavior  of  the  regions  generated  by  the  two  methods 
differs  with  the  number  of  stations  in  the  net,  configuration 
of  the  net,  standard  deviation  of  the  error  used,  and  relative 
location  of  target  with  respect  to  the  net.   Thus,  a  general 
case  may  not  exist.   A  valid  comparison  would  include 
sufficient  combinations  to  detect  all  possible  patterns  of 
behavior. 

A  useful  measure  of  effectiveness  for  the  regions  would 
be  the  difference  in  the  probability  level  in  each  of  the  two 
regions  occupied  by  the  true  position  of  the  target. 

The  standard  deviation  used  in  generating  bearing  error 
is  also  used  in  computing  BPE  and  confidence  regions.   A 
minor  modification  to  the  program  would  permit  the  examination 
of  the  consequences  of  incorrectly  estimating  the  actual 
bearing  distribution  of  a  station,  by  computing  these  with 
a  different  standard  deviation  that  that  used  to  generate 
the  error. 

The  errors  generated  are  from  a  normal  distribution,  and 
include  no  wild  bearing.   A  means  of  detection  for  wild 
bearings  is  not  included  in  the  program.   Such  a  routine 
may  be  added  to  the  program,  along  with  a  generation  routine 
for  wild  bearings. 
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The  modifications  above  would  permit  examination  of  more 
realistic  situations  through  the  use  of  operational  data, 
along  with  a  model  of  the  actual  net  which  supplied  the  data 
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V.   AN  ALTERNATE  APPROACH 

An  underlying  assumption  in  the  methods  discussed  is  that 
there  is  a  uniform  prior  probability  for  the  target  location. 
That  is,  the  mathematical  solution  of  the  BPE  and  confidence 
regions  are  inputs  to  the  solution  of  the  operational 
question  of  target  location.   Prior  knowledge  of  probable 
target  location  does  not  affect  the  solution.   This  approach 
is  intuitively  correct  generally. 

Operationally,  there  are  two  uses  of  position  finding 
information : 

(1)  To  determine  the  location  of  a  target  of  known 
or  unknown  identity; 

(2)  To  identify  a  target  from  a  set  of  known  targets 
by  its  location. 

This  second  use  introduces  prior  information  into  the 
problem,  and  occurs  usually  when  attempting  to  identify  a 
fixed  land-based  station. 

Butterly  [2]  describes  a  method  for  incorporation  of 
information  of  this  type  in  producing  a  Bayesian  position 
estimate,  though  he  implies  that  its  use  could  be  more  general. 

A  further  modification  of  the  program  would  permit  a  test 
of  the  efficacy  of  this  method  by  comparing  the  number  of 
correct  identifications  when  using  the  Bayesian  method  with 
the  number  when  selecting  the  target  in  the  lowest  level  of 
the  confidence  region.   Both  the  BVN  and  x   methods  could  be 
used  in  each. 
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THT  <;  FflRTRAN  PROGRAM  CALCULATES  RADIO  DIRECTION  FINDING  POSITION 
AND  CONFIDENCE  REGIONS.  SEVERAL  OPTIONS  IN  OUTPUT  ARE 
AND  ARE  DOCUMENTED  BY  COMMENT  CARDS  IN  THE  PROGRAM. 


ESTIMATES 
AVAILABLE, 


INPUTS  TO      g?gig^0ftRpiNDING  STATIONS  (ONE  CARD,  FORMAT  12) 
LONGITUDE,  AND  STANDARD  DEVIATION  OF  EACH  STATION 
EACH,  FORMAT  3F10.4) 
TARGETS  (ONE  CARD,  FORMAT  12) 
AND  LONGITUDE  OF  EACH  TARGET  (ONE  CARD  EACH, 


NUMBER  OF 
LATITUDE, 
(ONE  CARD 
NUMBER  OF 
LATITUDE 


FORMAT  2F10.+)  „  „_,_,. 

SIZE  OF  GRAPHICAL  PLOTS  IN  DEGREES 


(ONE  CARD,  FORMAT  F10.4) 


LIMITS 


-    STATIONS 
TARGETS 


8 


SIGN    CONVENTION:     LATITUDES 

LONGITUDES 
BEARI NGS 


NORTH  0  TO  +90.0 
SOUTH  0  TO  -9  0.0 
EAST  0    TO    -180.0 

WEST  0    TO    +180.0 

EAST    OF    NORTH    0    TO 
WEST    OF    NORTH    0    TO 


+180.0 
-180.0 


PART  ONE  -  INITIALIZATION  AND  INPUT 


GAMMA ( 8,8) 
STALON(IO) 
THETA(  10, 


DIMENSION  BETA(IO.IO) 
DIMENSION  STEEST (8, 8) 
DIMENSION  STALAT(U) 
DIMENSION  TARLON(  20) 
DIMENSION  WORM( 13,13) 
DIMENSION  PJ(1 0) ,ZETA(10) 
DIMENSION  Bl  (8) ,  B2(8)        „,-.„. 
DIMENSION  At (8) ,BE(8) ,CE(8),DE(8), 
DIMENSION  Y( 13,13) 
DIMENSION  S0RAD(8) 
DIMENSION  C(13,13) 
REAL*4  T(24) 


,  STDEV(IO),  TARLAK20) 
20)  ,  DISK  10,20) 


EE(8) ,FE(8) 


RADDEG=. 017453293 
IX=46721 
PI=3. 141592 
RADIUS=10803.0/PI 
EX    =    0.0 


READ    STATION    POINTS 
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READ(5, 100)  NSTA 
100  FGRMAT1I2) 

2  READ  (5*102)ASTALAT(I>,  STALCN(I),  STDEV(I) 
102  F0RMAT(3F10.4) 

READ  TARGET  POINTS 

READ(5»100)  NTAR 
DO  1  1=1 »NTAR 
1  READ(5,1'02)  TARLAT(I)»  TARLON(I) 

READ  GRAPH  SIZE 

READ(5,102)  CGRAPH 


PART  TWO  -  GENERATION  OF  BEARINGS  AND  BEARING  VARIANCES 


THIS  SECTION  OF  PART  TWO  COMPUTES  BEARING  ANGLES  AND  DISTANCES 
USING  NAPIER'S  ANALOGIES 

DO  10  1=1, NSTA 
DO  1C  J=1tNTAR 
AA=ABS(  STALGNU  J-TARLON(J)  ) 

IFUA.GT  .130.0  )  AA  =  360.0-AA 
GAMMA( I , J)=AA 

CHECK  FOR  SAME  MERIDIAN  CASE 

IF(AA.GT.O.Ol)  GO  TO  105 

DD-n    Q 

IF^STALAT(I)  .GT  .TARLAT(J))  BB= 180.0 

DIST  lit  J)  =  ABS(  STALATU  )-TARLAT(  J)  )  *RAD  I  U  S*RADDEG 
GO  TO  23 

105  CONTINUE 

CC  =  ABS( AA-180.  0) 
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c 
c 
c 


c 

c 
c 
c 


106 


C 

c 
c 
c 
c 
c 
c 


IFCCC.GT. 0.301)     GO    TO    106 

DD=0.0 

SUMLAT=STALAT(  I  J +T-ARLAT  {  J  ) 

IFCSUMLAT.LT. 0.0)  DD=180.0 

THETA( I ,  J)=DD 

EED=180.0-ABS(SUMLAT) 

DIST( I » J  )=EED*RADI US*RADDEG 

GO   TO    23 

CONTINUE 


CONTINUE    SOLUTION    OF    TRIANGLE 

WVALUE    =    GAMMMI  »J)*.5*RADDEG 
C0T1    =    COS(WVALUE)/SIN(WVALUE) 
BRAVC    =    93.0    -    TARLAT(J) 
ABLE    =    90.0    -    STALATU  ) 

C0S1  =  COS( (ABLE    -  BRAVO) *. 5*RADDEG ) 

C0S2  =  CGSUA3LE    +  BRAVC  )*. 5*RADDEG) 

SIN1  =  SIN((ABLE    -  BRAVO  )*.  5*RADDEG ) 

SIN2  =  SINUABLE    +  BRAVO)*. 5*RADDEG) 

YA=C0T1*C0S1 
YB=C0T1*SIN1 

ATA=ATAN2(YA,CCS2) 
ATB=ATAN2( YB ,SIN2) 
B3B=f ATA-ATBJ/RADDEG 


8BB=ANGLE,  NOW  COMPUTE  PROPER  SIGN 


SL=STALON( I ) 

TL=TARLON( J ) 

IF(SL.LT.O.O) 

IF(TL.LT.O.O) 

TL=TL-SL 

IF(TL.LT.O.O) 

IF(TL.LT. 180.0)  BBB=-BBB 

THETAt I , J)=BBB 


SL=360.0+SL 
TL=360.0+TL 

TL=TL+360.0 


THETA=TRUE  BEARING  FROM  STATION  TO  TARGET 

DISTANCE  FROM  STATION  TO  TARGET  USING  LAW  OF  COSINES 
GAMMA(ItJ)  =  GAMMAUt  J)*RADDEG 


17 


TARLAT( J)  =TARLAT(  J)*RADDEG 
^f|pAI(SIN(}ARLTAT(jn*sfN<isTALAT(I))+COSITARLAT(J))*COS(STALAT(I)) 

1*C0S( GAMMA (  I  J)  ) 
DFDIS'f     =    ARGOS  (TEMP) 
DIST(I,J)     =    DFDIST*RADIUS 

DIST=DISTANCE    IN    NAUTICAL    MILES 

TARLATC J)=TARLAT( J) /RADDEG 
STALAT(I)=STALAT(I)/RADDEG 

23    CONTINUE 


C 

C 


THIS    SFCTION    OF    PART    TWO    ADJUSTS    THE    STANDARD    DEVIATION    TO    BE    USED 
!m    r-Jrtnii-     FPQHR    COMPUTATION    POR    THE    DISTANCE     BETWEEN    STATION    AND 
TARGET^       THE    SOURCE "oFTHI S    COMPUTATION    IS    NOT    KNOWN. 

SD    =    STDEV( I ) 

R    =    DIST( I, J )/ 100.0 

IF    (R-54.0)     20,20,22 

22  SD    =    SD*     (R/(108.0    -    R) ) 
GO    TO    3  2 

23  IF     (R-1J.0)     2^,26,26 

24  IF     (R-4.C)     28,30,30 

28    SD    =    PI/2.0 

GO    TO    32  ,  r  ^    n  n    m 

30    SD    =    SD*(.0204*R*R    -    .4C2*R    +    3.0) 

26    SD    IsD*(.0a071*R*R    -    .3213*R    +    1.1598) 
32    STFESTiUJ)    =    SD*SD*RADDEG*RADDEG 

10    CONTINUE 

253    F0RMAT6(///!36X,45H    BEARING    ANGLE    FROM     ITH    STATION    TO    JTH    TARGET,/ 

1WRITE<6,256)((I,J,THETA(I,J),J=1.NTAR) ,I=1,NSTA) 
256    F0RMAT(39X,13HFRGM    STATION     ,12, 
111H    TO    TARGET     ,I2,4H    —     ,F10.4J 

255    F0RMAT6(/A40X,40HDISTANCES    FROM    ITH    STATION    TO    JTH    TARGET,//) 
WRITEC6,256)(II,J,DIST(I,J),J=1,NTAR),I=1,NSTA) 


THIS  SECTION  OF  PART  TWO  COMPUTES  EEARING  ERROR  FROM  A  SIMULATED 
NORMAL  DISTRIBUTION  USING  A  RANDOM  NUMBER  GENERATOR , AND  ADDS  THE 
ERROR  TO  THE  TRUE  BEARING, 

WRITE(6,260) 
260  FORMAT  { // ,50X , 1 5HNGRMAL  VARIATES,//) 

DO    40     I=1,NSTA 

DO    40    J=1,NTAR 

STEESTU.Jl     =    STEESTU  » Ji/RADDEG 

CALL  GASS( IX,STEEST(I , J) ,EX,X) 

BETAU  ,  J)=THETA(  I ,  J ) +X 

WRITE(6,#256)I,  J,X 

40  CONTINUE 

OUTPUT  ANGLES  WITH  ASSOCIATED  ERROR 

WRITE(6t200) 
200  FORMAT  ( //6X , 6H STA LAT ,  8X,  6HSTAL0N,  8X,6HTARLAT,  8X,6HTARL0N, 
18X,10HTHETAU,  J  )  ,  4X,  5HSTDEV) 

WRITE (6,259) ({ STALATII ) ,STALON( I ) »TARLAT (J) ,TARLON( J) , 
1  BETA( I , J) ,STEEST{ I ,J) ,  J=1,NTAR),  I=1,NSTA) 
259  FORMAT {/,6 (F12.5.2X )) 


PART    THRE 


COMPUTATION    OF     POSITION    ESTIMATE,     OR    FIX 


DO    3321     J1=1,NTAR 

WRITE(6,220)    TARLAT(Jl),     TARLON(Jl) 

220  F0RMAT(71H1  THE    FOLLOWING    BEARINGS    WERE    TAKEN    ON    A    TARGET     LOCA 
1TED    AT    LATITUDE,     F12.5,11H    LUNG  I  TUDE , F 1 2. 5 , // ) 

WRITE(6,300) 
300    FORMAT ( 10X, 8HL ATI TUDE  ,  12 X , 9H LONG  I TUDE , 14X , 5HTHETA , 1 IX , 1 2 HTRU E     BEAR 
1ING,     10X,     7HSTD    DEV,/J 

WR[TE(6,221)     (     STALATU  )  ,STALON{  I  )  ,     BET  A  (  I ,  Jl  )  ,  THET  A  (  I  ,  J  1  ) , 
1STDEV(  I )     ,  1  =  1, NSTA) 

221  F0RMAT(5(10X,F10.5) ) 
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THIS    SECT  ION    OF    PART 
THE    FIX    USING    POPE'S 


THREE    CALCULATES    A    FIRST 
VECTOR    METHOD    (1971 ) 


POINT    ESTIMATE    OF 


COMPUTE    BEARING    PLANE     INTERSECTIONS 

WRITE(6,701) 
701     FORMAT { // ,44X, 26HBEARI NG    LINE    INTERSECTIONS, 
1/,40X,8HSTATI0NS»4X,8HLATITUDE,4X. 
29HL0NGITUDE,// ) 
DO    88    I=1,NSTA 


Bid)    = 
B2(I)     = 

SLAT 
CLAT 
SLON 
CLON 
SBNG 
CBNG 


STALAT( I ) 
STALONl I ) 

SINCBK  I  )  ) 
COS (Bl(  I)  ) 

■SIN(B2(  I  )  ) 
C0S(B2(  1)  ) 
SIN(BETA(  I 
COS(BETA(  I 


*  RADDEG 

*  RADDEG 


J1)*RADDEG) 
JD-'RADDEG) 


XE  =  CLAT  *  CLON 
YE  =  CLAT  *  SLON 
ZE    =    SLAT 

AF(I)    =    CBNG*SLON    -    SBNG* SLAT* CLON 
BEtl)     =    -CBNG*CLON    -    S BNG*SL AT*SLGN 
CS(I>     =    SBNG-CLAT 


DE(I)  =  BE( I)*ZE 
EE(I )  =  CE( I)*XE 
FECI)    =    AE(I)*YE 

.3    CONTINUE 

USUM  =  0.0 
VSUM  =  3.0 
WSUM    =    0.0 

DO    89    I=2,NSTA 

iLiniT  =   i-i 

DO    90    J=l, ILIMIT 

UE  =  3E( I )*CE( J  ) 
VE  =  CE( I J*AE( J) 
WE    =    AE( I)*BE( J) 


CE( I )*YE 
AE( I)*ZE 

BE( I )*XE 


CE( I)*BE( J) 
AE( I)*CE(J) 
BE(I)*AE(J) 
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TI  =  UE*DE(I)  +  VE*EE(I)  +  WE*FE(I) 
TJ  =  UE*DE(J)  +  VE*EE(J)  +  WE*FE(J) 

IF  (TI*TJ  .LT.  0.0)  GO  TO  90 
IF  (TI  .GT.  0.0)  GO  TO  91 

UE  =  -UE 
VE  =  -VE 
WE  =  -WE 

COMPUTE  CENTROID  OF  INTERSECTIONS 

91  USUM  =  USUM  +  UE 
VSUM  =  VSUM  +  VE 
WSUM  =  WSUM  +  WE 

UMAG=SQRT(UE*UE+VE*VE+WE*WE) 

ULAT=ARSIN(W5/UMAG) /RADDEG 
UL0N=-ATAN2(VE,UE) /RADDEG 

WRITE  (  6,702)  I  ,  J,ULAT,ULCN 
702    F0RMAT(40Xt I2t5H    AND    , I  2 , 3X, F 8. 3 , 4X, F8. 3 ) 

ULAT,     ULON    =    BEARING    PLANE     INTERSECTIONS 

90    CONTINUE 
89    CONTINUE 

CONVERT    CENTROID    TO    LATITUDE    AND    LONGITUDE    OF     FIRST     POINT     ESTIMATE 

X1L0N    =-ATAN2( VSUMtUSUMJ/RADDEG  „ 

X1LAT=ARSIN( wSUM/SQRT ( USUM*USUM+VSUM* VSUM+WSUM^WSUM ) ) /RADDEG 

WRITE(6f222) 

WRITE (6,223)     X1LAT,     X1L0N 

222  FORMAT!//, 10X.24HTHE    FIRST    POINT    ESTIMATE,//) 

223  F0RMAT(10X,8HX1LAT    =     , F10.5 t 5X, 8HX1L0N    =     ,F13.5t//) 

THIS    SECTION    OF    PART    THREE    TRANSFERS    THE    PROBLEM    TO    A    PLANE    TANGENT 
TO    THE    EARTH    AT    THE    FIRST    POINT    ESTIMATE,     THEN    COMPUTES    A    POSITION 
ESTIMATE    USING    THE    LEAST    SQUARES    METHOD    OF    KUKES    AND    STARIK. 

ASUM=0.0 
BSUM=0.0 
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CSUM=0.0 
DSUM=Q. 0 

ESI'M  =  0.0 

DO  8999  L=1.NSTA 
SDRAD(L)=STDE\ML)*RADDEG 
8999  CONTINUE 

DO  900  2  I=1,NSTA 
XSI=BETA( I, J1)*RADDEG 
XSI=ABS{XSI ) 

XSI=STATION  BEARING  IN  RADIANS 

COMPUTE  BFARING  IN  LOCAL  REFERENCE  PLANE  USING  NAPIER'S  ANALOGIES 
AND  LAW  OF  SINES 


RHO=(STALON( IJ-X1L0N 
RH0=A3S(RH0) 

IF1RHG.GT.PI )     RH0=(2.0*PI) 


}*RADDEG 

•RHO 

CHECK    FOR    DUE    NORTH    OR    SOUTH    BEARING 

IFtRHO.GT. 0.001)     GO    TO    345 

IF(STALAT(  D.LT.X1LAT)     ZETA(I) =0.0 
GO    TO    347 

345  CONTINUE 
RR=ABS(RHO-PI) 
IFtRR.GT. 0.001)     GO    TO    346 
ZETACI)=PI 

P0LE=STALAT(I)+X1LAT 
IFIPOLF.LT.0.0 )     ZETA(I)=0.0 
GO    TO    347 

346  CONTINUE 

CONTINUE     SOLUTION     OF    TRIANGLE 

CEE=(90.-STALAT(I) )*RADDEG 
DENOM=.  5*(  XSH-RHO) 
RUM11=.5*(XSI-RH0J 

RUM12=. 5*CEE 

ANUM=SIN(RUM11 )*TAN(RUM12) 
DENA=SIN(DENOM) 
ARKTN1=ATAN2(ANUM,DENA) 
RTA=ABS(  DENOM-PI/2) 
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RTB=ABS(RUMll-PI/2)  _ 
IFIRTA.GT. 0.001)  GO  TO  355 
BNUM=C0S(RUM11)*TAN(RUM12) 

DENB=0.001 

355  IFtRTB.GT. 0.001)  GO  TO  356 
BNUM=0.001 

DENB=COS(DENOM) 
GO  TO  357 

356  BNUM=C0S(RUM11)*TAN(RUM12) 
DENB-COS(DENGM) 

357  CONTINUE 
ARKTN2=ATAN2(BNUMtDENBJ 
AAEE=ARKTN1+ARKTN2 
IFIXSI.LT. 0.001)  XSI=0.001 
CDD=SIN(CEE)*SIN(XSI)     ^^n 
1F(AAEE.LT. 0.001)  AAEE=0.001 
EEE=SIN( AAEE) 

FFF=DDD/EEE  ,™„«„ 

IF(FFF.GE.l.O)  FFF  =  0. 999999 
CCEE=ARSIN(FFFJ 

C 

ZETA(I)=PI-CCEE  ^^^^ 

IF(AAEE.LT.CEE)     ZETAU)=CCEE  ,c_.fT1 

IF(BETA{ I, Jl) .LT.O.O)     ZETA ( I ) =-ZETA ( I ) 
C 

347    CGNTINUE 
C 
C       ZETA=L0CALIZED    BEARING 

C 

C   TAKE  NEGATIVE  OF  ZETA  TO  CONVERT  TO  NORMAL  TRIG  SIGN  CONVENTION 

C 

ZETAU  )  =-ZETA(  I) 
C 

C   CALCULATE  P FRPFMD ICULAR  DISTANCE  FROM  LOCALIZED  BEARING  TO 

C   CENTER  OF  COORDINATE  SYSTEM 

C 

YPRIM=(90.-X1LAT      )*RADDtG 

AAEE=ABS(AAEE) 
YDIST=    (YPRIM-AAEF) 
PJ(  I )=YDIST^SIN(ZETA(I  )) 
C 

C   SUM  TERMS  FOR  FOLLOWING  CALCULATIONS 
C 
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ASUM=ASUM+(COS(ZETA(I ) ) *COS ( ZET A { I ) ) ) / (SDRAD1 I )*S CRADj  I  )  ) 
BSUM=BSUM+    SIM     ZETA(I))*COS(ZETAiI}))/(SDRAD(I    *SDRAD  (I 
CSUM=CSUM+1SIN(ZETA(I))^SIN{ZETA{I)})/(SDRAD(I)*SDRAD(I)) 

C 

9002  CONTINUE 
C 

DSW=DSUM+PJI(,kI*{(BSUM*SIN(ZETA{K))-CSUM*COS(ZETA(K)))/(SDRAD(K)* 

1ESUM£ESUMipJ(K)*((ASUM*SIN(ZETA(KJ)-BSUM*COS(ZETA(KJ ) ) / (S DR AD { K )* 
1SDRAD(K) 3 ) 
C 

9006  CONTINUE 
C 

C      CALCULATE    LEAST    SQUARES    ESTIMATES    AND    CONVERT    TO    DEGREES 

C0EF=1.0/( (ASUM*CSUM)-(BSUM*BSUM) ) 

XO=COEF*DSUM 

YO^COEF*ESUM 

S0KO  =  -2.*AL0G(  .1  ) 

XEST=X1L0N+X0/RADDEG 

YEST=X1LAT+Y0/RADDEG 

C       XEST,    YEST    =    LEAST    SQUARES    ESTIMATE    OF    POSITION 

C 
C 

9007  F0RMAlt'9??^THEELEAST  SQUARES  ESTIMATE  OF  POSITION  IS  •  ,F6. 2 , 
1«    DEGREES  LATITUDE'  //43X,F7.2,»   DEGREES  LONGITUDE '/// ) 

C 

c 
c 
c 

r 

C   PART  FOUR  -  COMPUTATION  OF  CONFIDENCE  REGIONS 
C 

c 

C   THIS  SECTION  OF  PART  FOUR  CALCULATES  A  BIVARIATE  NORMAL  CONFIDENCE 

C   REGION  (DANIELS,  KUKES  AND  STAR1KJ 

C 

C   CALCULATE  ANGLE  BETWEEN  MAJOR  AXIS  OF  ELLIPSE  AND  MERIDIAN 


C 


AA1=2.*BSUM 
AA2=CSUM-ASUM 

TW0ANG=ATAN2(AA1»AA2) 
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BANG=TW0ANG/2. 0 
C 

IF(BANG)651,651,652 

651  BANG=BANG+(PI/2.0) 

GO   TO    653 

652  BANG=BANG-(PI/2.0) 

653  CONTINUF 
ANG=BANG/RADDEG 

C 
C 

c 

C   CALCULATE  SEMIAXES  FOR  ELLIPSE  OF  .9  PROBABILITY 
C 

IF(BANG.GT.1.56J  BANG=1.56 

IF(BANG.LT.-1.56)  BANG=-1.56 
C 

SDENOM=CSUM+BSUM*TANlBANG) 

FALSE=ABS(SDENOM) 

SAMAJO=SORT( SOKO/FALSE) 

SDENO=ASUM-BSUM*TAN(BANG) 

FATS=ABS(SDENO) 
C 

SAMINO=SQRT(SQKO/FATS) 

S A  M I NO  =  S AM I  NO/ R A  CD  EG 

SAMAJO=SAMAJO/RADDEG 
C 
C 

c 

WRITE( 6, 9081)   ANG , SAMAJO , SAMI NO 
9081  FORMAT {«  »,10X,1 ANGLE  BETWEEN  MAJOR  AXIS  AND  MERIDIAN  =  «,F6.2// 
134X, 'SEMI-MAJOR  AXIS  =  •  , F6 .2 » // t 34X * ' S EM  1 -M I NuR  AXIS  =  ', 
C 
C 

c 
c 

C   THIS  SECTION  OF  PART  FOUR  PRINTS  A  GRAPHIC  PLOT  OF  A  BIVARIATE 

C   NORMAL  CONFIDENCE  REGION.   THE  SIZE  IS  DETERMINED  BY  THE  INPUT 

C   •CGRAPH*.   USE  THE  SMALLEST  VALUE  OF  CGRAPH  THAT  WILL  DISPLAY 

C   THE  REGION*  TO  MINIMIZE  DISTORTION. 

C 

C       THIS    SECTION    MAY    BE    SKIPPED    IF    THIS    OUTPUT     IS    NOT    DESIRED.        JUMP 

C       FROM    THIS    POINT    TO    THE    BEGINNING    OF    THE    CHI-SQUARE     SECTION,     BELOh 

C 

C       LOAD    PROBABILITY    ARRAY    FOR    BIVARIATE    NORMAL 

C 

XOR=XEST 

YOR=YEST 

DO    9311    J=l, 13 
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90 
90 


12 
11 


EXVAL=(CGRAPH/2.0)-(( ( J-l) /12. 0 ) *CGRA PH ) 

EXVAL=EXVAL*RADDEG 

DO    9012    1=1,13 

YVAL={CGRAPH/2.0)-(  ( ( 1-1 J / 12. 0 ) *CGRAPH) 

YVAL=YVAL*RADDEG 

H00=BSUM#EXVAL*YVAL-(ASUM/2. ) * EXVAL* EXVAL- (CSUM/2 .0 )*YV AU 

PROB=1.0-EXP(HOO) 

WORMU  , J)=PROB 

CONTINUE 

CGNT INUE 


:YVAL 
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N=13 
M=13 

AZ=0.0 

BZ=0.0 

BND=0.0 

AMIN=0.0 

IJT=0.0 

IC0N=1 

IGR  =  1 

DO  9080 

T(LL)=0 


LL=1,24 

0 


PRINT  THE  CONFIDENCE  REGION  PLOT 

CALL  MTRMAP ( WORM ♦ N, M, T,BND,AZ,BZ, AMI N,IJT,  ICON, IGR) 
WRITE(6,591) 
591  FORMAT(///,45X, ■ BIVARIATE  NORMAL  CONFIDENCE  REGION') 

WRITE (6, 509)  CGRAPH,CGRAPH 
509  FORMAT(///,20X,37H   THE  SIZE  OF  THE  CONFIDENCE  PLOT  IS  ,F10.4, 
114H  DEG.  LAT.  BY  .F1J.4,  12H  DEG.  LCNG.  ,/,40X, 
241H   THE  FIX  POINT  IS  LOCATED  AT  (007,  007).  ) 


THIS  SECTION  OF  PART  FOUR  CALCULATES  A  CHI-SQUARE  CONFIDENCE 
REGION,  AND  PRINTS  A  GRAPHICAL  PLOT.   NO  OTHER  OUTPUT  IS  AVAILABLE. 
ADDITIONALLY,  LEVELS  OF  CONSTANT  VALUE  OF  Q  MAY  BE  PLOTTED  (DANIELS) 
THIS  SECTION  MAY  BE  SKIPPED  BY  JUMPING  FROM  THIS  POINT  TO  THE  END  OF 
THE  MAIN  PROGRAM.    EITHER  THE  0  OR  CHI-SQUARE  PLOTS  MAY  BE  SKIPPED 
SEPARATELY  BY  JUMPING  AROUND  THE  APPROPRIATE  'CALL  MTRMAP...'  CARD 
AT  THE  END  OF  THIS  SECTION. 


X1LAT=YEST 
X1L0N=XEST 
X1LAT  =  X1LAT 
X1L0N  =  X1L0N 


CGRAPH/2.0 
CGRAPH/2.0 
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610 


620 


571 
502 
501 


DO  501  K  =  1,13 
DO  502  L  =  1,13 
Y(K,L)    =    0.0 


DO    571 

G=ABS( S 

IFtG.GT 

IFtG.GT 

ANGLE=0 

IF(STAL 

GO    TO    6 

CONTINU 

ABLE=90 

BRAV0=9 

WVALUE= 

C0T1 

C0S1 

C0S2 

SI  Ml    = 

SIN2    = 

RR=C0T1 

SS=ATAN 

TT=C0T1 

UU=ATAN 

ANGLE=( 

CONTINU 

SL=STAL 

TL=X1LC 

IFtSL.L 

IFtTL.L 

TL=TL-S 

IF(TL.L 

IFtTL.L 

CHIDIF= 

Y(K,L) 

PAR1=NS 

PAR2=Y( 

CALL    GA 

C(K,L)= 

CONTINU 

X1LCN  = 
X1LCN  = 
X1LAT    = 


1=1, 

TALO 

.180 

.0.0 

.0 

AT(  I 

20 

E 

.0-S 

0.0- 

G*.5 

COS( 

cost 
cost 

S I  N  ( 
SI  N{ 
*SIN 

21RR 
-COS 
2(TT 

uu-s 


NSTA 

N( I J-X1L0N) 
.0)     G=360.0-G 
01)     GO    TO   610 

J.GT.X1LAT)     ANGLE=180.0 


TAL 

X1L 

*RA 

WVA 

(A3 

(AB 

(AB 

(AB 

1 

,SI 

1 

,CD 

S)/ 


AT  (I  J 

AT 

DDEG 

LUE)/SIN(WVALUE) 

LE-BRAVO)*.5*P.ADDEG) 

LE+BRAVD)*.5*RADDEGJ 

LE-BRAVO J*.5*RADDEG) 

LE-t-BRAV0)*.5*RADDEG) 

N2) 

S2) 

RADDEG 


CN( 

N 

T.O 

T.O 

L 

T.O 

T.l 

ANG 

=  Y 

TA/ 

K,L 

MA 

1.- 

E 


I) 

.0)  SL=360.0+SL 
.0)  TL=360„0+TL 

.0)    TL=TL+360.0 

t.3.0)     ANGLE=-ANGLE 

LE-3ETM  I  ,J1  ) 

(K,L)  +  CHIDIF*CHIDIF/STDEV(I }/STDEV( I ) 

2. 

)/2. 

(PARI, PAR2, GAM, B, 0.0) 

GAM 


X1LCN    -    CGRAPH/12.0 
X1LGN    +    CGRAPH* 13. 0/12.0 
X1LAT    -    CGRAPH/12.0 


SET    UP    THE    CONTOUR    SUBROUTINE 
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N  =  13 

M  =  13 

AZ  =  0.0 

BZ  =  3.0 

BND  =  0.0 

AMIN  =  0.0 

IJT  =  0 

ICON  =  1 

IGR  =  1 

DO  504  II  =  1,24 

504  T(  I  I  )  =  0.0 

c 

, 

c 

PLOT  THE  0  LEVELS 

c 

CALL    MTRMAPtY, N, M, T , BND , AZ , BZ , AMIN, IJT, ICON  ,  IGR ) 
WRITE{6,592) 

592  F0RMAT(///,45X, » LEVELS    OF    CONSTANT    VALUE     OF     Q» ) 
WRITEi 6,509)     CGRAPH  ,CGRAPH 

PLOT    THE    CHI-SQUARE    REGION 

CALL    MT RMAP ( C,  N, M, T, BND, AZ,BZ, AMIN, IJT, ICON, IGR) 

WRITE(6,593) 

593  FORMAT! ///,45X, 'CHI-SQUARE  CONFIDENCE  REGION') 
WRITE(6,f09)  CGRAPh,  CGRAPH 

3321  CONTINUE 

STOP 
END 


SUBROUTINE    GASS ( I  X , S, AM, V) 
A    =    0.0 
DO    50    1=1.12 
CALL    WRANDU( IX, IY, Y) 
IX    =    IY 
50    A    =    A+Y 

V    =     (A-6.0)=-S  +  AM 

RETURN 

END 
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SUBROUTINE     WRANDUt I  X, I Y »Y) 
IY    =     IX*65539 
IF     (IY)     5,6,6 

5  IY    =    IY    +    2147483647+1 

6  Y    =     IY 

Y    =    Y*.4656613E-9 

RETURN 

END 
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